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A self-homodyne detection scheme is proposed to perform 
two-mode tomography on a twin-beam state at the output of 
a nondegenerate optical parametric amplifier. This scheme 
has been devised to improve the matching between the local 
oscillator and the signal modes, which is the main limitation 
to the overall quantum efficiency in conventional homodyning. 
The feasibility of the measurement is analyzed on the basis 
of Monte-Carlo simulations, studying the effect of non-unit 
quantum efficiency on detection of the correlation and the 
total photon-number oscillations of the twin-beam state. 



I. INTRODUCTION 

One of the most significant advances in modern quan- 
tum optics is the theoretical development Gl and subse- 
quent experimental realization Q of homodyne tomog- 
raphy. This measurement scheme allows one to recon- 
struct the density matrix of the quantum state from a set 
of field quadratures measured by a balanced homodyne 
detector. Reconstruction methods, initially based on ap- 
proximate inverse Radon transform of the quadratures 
histograms, have been enhanced later through exact algo- 
rithms |3|-|6|] that achieve the measurement of the matrix 
element by sampling a corresponding pattern function of 
the experimental homodyne outcomes (for a review see 
0). These algorithms have been proven to be very sta- 
ble and fast enough to allow real-time data sampling. 
For the photon-number representation, the calculation 
of the pattern functions has been greatly improved by 
means of factorization formulas || and asymptotic ap- 
proximations H for large photon numbers of the matrix 
indices. The direct sampling approach has been imple- 
mented experimentally to measure the photon statistics 
of a semiconductor laser []l0|| , and the density matrix 
of a squeezed vacuum pT| . The success of optical ho- 
modyne tomography has stimulated the development of 
state- reconstruction procedures for atomic beams [ ^2| , 
the experimental determination of the vibrational state 
of a molecule Jl3) , of an ensemble of helium atoms juj , 
and of a single ion in a Paul trap Q. Finally, some 
non-tomographic state reconstruction methods have also 
been recently proposed p^ ]. 

While the full density matrix reconstruction requires 
the knowledge of the phase of the detected mode with 
respect to the local oscillator (LO), for the diagonal 



matrix elements it is just sufficient to average over a 
random phase fTc}] . The typical nonclassical states of 
interest — such as squeezed states — already exhibit inter- 
esting quantum features in just the photon number dis- 
tribution; this makes homodyne tomography especially 
attractive. Among the quantum features of interest, 
there are the even-odd oscillations in the photon num- 
ber distribution of a squeezed vacuum flr7| , which were 
recently observed experimentally |Tl[| . In two-mode to- 
mography of a twin-beam state produced in parametric 
down-conversion, we are interested in features of the joint 
photon-number distribution of the signal and the idler, 
such as the delta-function correlation between the photon 
numbers of the two modes, and the even-odd oscillations 
of the total photon number. The sampling algorithm 
for the two-mode tomography is obtained by a straight- 
forward extension of the single- mode case |l8|j . In the 
relatively new field of multimode tomography, recent ad- 
vances have been made in the theoretical description fis| ] 
and the experimental measurement [ppf of the photon- 
number correlation between two temporal modes. 

From the experimental point of view, homodyne to- 
mography of the photon-number distribution is a viable 
alternative to direct detection. It allows one to measure 
very weak photon fluxes — of the order of a fraction of a 
photon per measurement time — using high quantum ef- 
ficiency fast p-i-n photodiodes, as compared to the slow 
and less efficient avalanche photodiodes used for direct 
detection. This convenience, however, comes with its own 
price tag. One encounters the problem of mode matching 
between the LO and the detected modes ^lj , determined 
by their spatial/temporal overlap, which gives a detri- 
mental contribution to the overall quantum efficiency. As 
shown in Ref. jij , the detection of the quantum features 
is rapidly degraded by less-than-unity quantum efficiency 
of the homodyne detector, and the degree of degradation 
rapidly increases for larger photon numbers. 

The problem of mode matching becomes especially 
severe for quantum states generated in traveling-wave 
or pulsed experiments, particularly those employing the 
optical-parametric amplifiers (OPA's). It has been shown 
that a LO well matched to a squeezed vacuum can be gen- 
erated in the same parametric process 22 . For example, 
in Ref. [^2| a polarizationally nondegenerate OPA was 
used to produce the squeezed vacuum and the matched 
LO in two orthogonal polarizations. However, while this 
approach is justified for measurements of the squeezing, 
it cannot be used for density-matrix reconstruction, as a 
tiny leakage of light from the LO polarization into the 
signal polarization can easily spoil the signal photon- 
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number distribution. 

In this paper, we address the problem of generating a 
matched LO for the reconstruction of the density matrix 
of the output state of a polarization-and-frequency non- 
degenerate OPA. In the spirit of Ref. p2[ , we develop a 
concept of self-homodyning that allows one to create both 
the LO and the signal in the same OPA. In the direct de- 
tection of the output signal field, a strong mean field at 
the central frequency u>o can serve as a LO for measuring 
a mode that consists of two sidebands at ujq ± O. As we 
will show in the following, the relative phase between the 
LO and the two-sideband mode can be varied, thus al- 
lowing homodyne tomography of the latter. In this way 
one can perform the tomographic reconstruction of full 
joint density matrix of the signal and idler twin-beam 
modes. In this paper, we consider the measurement of 
the joint photon-number distribution of these two modes 
and the photon-number distribution of the signal mode 
alone. For the latter, a thermal distribution is expected, 
as seen in recent self-homodyning experiments ]23| ]. We 
also show that self-homodyning can be used to measure 
the photon statistics of the +45°- and the — 45°-polarized 
linear combinations of the signal and idler modes. 

From Monte-Carlo simulations we will estimate the 
experimental conditions that are needed to extract the 
joint photon-number probability distribution of the twin 
beams, the photon correlation between the modes, and 
the quantum oscillations of the total photon number. 
We will show how these quantities can be experimen- 
tally measured for realistic values of quantum efficiency 
(^0.9) of the photodiodes and for reasonable number of 
data points (~ 10 6 ). 

In Section || we give a detailed theoretical description 
of the self-homodyne measurement, relating the measure- 
ment of the fie ld qu adratures to the output photocurrents 
in Subsection HA, and evaluating the joint pro babi lity 
distribution of the photocurrents in Subsection II B, in 



tocurrent distribution used in Subsection II B 



a form suitable for Monte-Carlo simulations, also taking 
into account the effect of non-unit quantum efficiency. In 
Section [II we briefly review the exact reconstruction al- 



gorithm for quantum tomography, for one mode only in 
Subsection III A, and th en wi th extension to any n umber 
of modes in Subsection [II B . In Subsection III Cj we an- 
alyze how the two-mode tomography is a chieve d through 
self-homodyne detection. In Subsection III D| we intro- 
duce the concepts of the measurement of the "dressed" 
state, often adopted in experiments, — as opposed to the 
"bare" state, usually assumed by the theorists. In Sec- 
tion |^ we present some selected Monte-Carlo simula- 
tions, also for non-unit quantum efficiency, for both the 
bare and the dressed states. We will focus attention on 
the joint photon-number probability, on the correlation 
between the photon numbers of the two modes, and fi- 
nally on the total photon-number probability, which ex- 
hibits oscillations typical of the twin-beam state. Section 
^ concludes the paper with a discussion of the results in 
view of the feasibility of the real experiment. The Ap- 
pendix covers the details of derivation of the joint pho- 



II. 



THEORETICAL DESCRIPTION OF THE 
SELF-HOMODYNE MEASUREMENT 

A. The detector 



The scheme of a self-homodyne detector is depicted in 
Fig. [I], along with the relevant modes of the electromag- 
netic field involved in the measurement. A nondegenerate 
optical parametric amplifier (NOPA) is injected with an 
input field having a strong coherent component at fre- 
quency luo with amplitudes aj and depending on the 
polarization, | denoting the vertical and <-» the horizontal 
polarization, respectively. The amplifier is pumped at the 
second harmonic cu p = 2ujq with amplitude a p 3> aj , a H , 
such that the pump can be considered as classical and un- 
depleted during the amplification process. At the output 
of the amplifier two photodetectors separately measure 
the intensities of a couple of orthogonally-polarized com- 
ponents of the field £j and £^. At the output of the 
photodetectors, a narrow band of the photocurrent is se- 
lected, centered around frequency f2 <C loq (typically ojq is 
optical/infrared, whereas fHs a radio frequency). In the 
narrowband approximation and for radiation absorbed in 
a thin detector layer, the filtered output photocurrents 
are given by the (complex) operators 
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dwf-(w + n)f+H, 7T = {!,<-}, (i) 

where :: denote the customary normal ordering with the 
(output) field annihilation-operator components £^ on 
the right and the creation operators £ ~ on the left, and 
the subindex 7r runs on the two independent polariza- 
tions | and <->. In terms of the annihilation and creation 
operators 8 and 8^ of the relevant output modes one has 



4(0) = 8J W S_* + 8^80* , 



(2) 



where the subindex refers to the central mode at fre- 
quency ujq, and ± refer to the sidebands at frequencies 
lu ± O, respectively. The input-output Heisenberg evo- 
lutions of the relevant field modes across the NOPA are 
given by 

80I = M«ol + , 8 ^ = ^a ^ + ^ ajj , 

8+j = /ia+| + va}_^ , 8 + ^ = /ia + ^ + va}_^ , (3) 

8_j = /ia-j + , 8_^ = /ia_^ + t^a^i , 

where a and a* denote the annihilation and creation op- 
erators for the input modes, /i = coshr, v = e l6p sinh?*, 
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FIG. 1. Scheme of a self-homodyne detector along with 
the relevant modes of the electromagnetic field involved in 
the measurement. The nondegenerate optical parametric am- 
plifier (NOPA) is seeded with input fields having strong co- 
herent components at frequency u>o, and is pumped at the 
second harmonic uj p — 2uo. At the output of the amplifier 
the intensities of the two different polarizations are separately 
measured by photodetectors, and a narrow band of the output 
photocurrents is selected, centered around frequency £1 -C u>o- 

r oc Lx^lo^l (L is the amplifier length, is the effec- 
tive second-order susceptibility). In the following we put 
8 p = 0, namely we set the pump phase as the reference 
phase for all modes. We assume the mode Sot to be in 
a highly excited coherent state with amplitude evj. For 
the purpose of measurement of the joint photon-number 
distribution, the mode clq^ will also be assumed in a 
highly excited coherent state with amplitude . In the 
case where we are interested in measuring the photon- 
number distribution of one beam only, the photocurrent 
produced by the '<-> '-polarized beam can be ignored or 
the mode can be assumed to be in the vacuum state. 
In the process of direct detection, the highly-excited cen- 
tral modes <2ot,<-> beat with the ujq ± ft sideband modes, 
thus playing the role of the LO of homodyne and hetero- 
dyne detectors. This converts the direct detectors into 
self-homodyne detectors whose experimental outcomes 
are the measured values of the following rescaled output 
photocurrents in the limit of strong LO's: 



ij(n) 



Tr L o[/T(fi)/5 L o] 
lim —= 
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Tr LO [/«(fi)/5 L o] 
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(4) 
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where ryj and denote the quantum efficiencies of the 
two photodetectors, a* denotes the complex conjugate of 
a, plo represents the density operator of the LO state, 
and Ttlo denotes the partial trace over the LO modes. 
In Eq. (||) 7 w (f2) is modified from that in Eq. (||) because 
of the non-unity quantum efficiencies of the two photode- 
tectors. It is given by 



(5) 



where b' = y/rj^ b 17I + ^/l — r)^ v 17I . Here v 17I for 7 S 
{0, +, — } and 7r 6 {|,^} are independent vacuum-state 



operators accounting for the loss at the three frequency 
components of each polarization mode. For the sake of 
simplicity, we will assume rjt = r}<^ = 1 for the rest of this 
Subsection. We will take into account the effect of non- 
unity quantum efficiency o n th e photocurrent probability 

Thus, using Eqs. 



If (Q) distribution in Subsection [I B 



3|) and 



quantum efficiency equal to unity, one obtains 



i:(0) = i(e-^_ I + e^ I 
i«(0) = -^( e -^S_ M + e^ 



(6) 



where <f> = arg(evj + ra^) is the phase of the mode 6 j 
relative to that of the pump with t = tanhr, and analo- 
gously tp = arg(a<_> + Ta j)- Taking the real part of the 
photocurrents at given radio-frequency phases £ and \ 
one has 



Re 



«1 = 



(7) 



where the operator X ,/,(£) denotes the quadrature at 
phase (f> of the mode with annihilation operator c, namely, 



(8) 



and the operator B\ is the annihilator of the polarized 
output mode 



(9) 



It is easy to check that the output modes (|^) have corre- 
sponding input modes given by 



=^=(e iX a^ + e- iX a + ^) , 
and they are related by the Heisenberg evolutions 



(10) 



(11) 



By scanning the relative phase cj> between & j and the 
pump mode, one can measure any quadrature X^[B^ ) 
of the output field. If the input sideband modes a± are in 
a state with a completely random phase, such as the vac- 
uum, then the only phase reference in the output modes 
b± is the pump phase 8 p — 0. In that case, the phase 
4> can be easily changed by delaying all the input fields 
with respect to the pump field, with no need to change 
the phase of Sot separately from the other input modes. 

From Eq. (||) one can recognize that there are actually 
four output modes that commute with each other; hence, 
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their quadratures could be jointly measured by the self- 
homodyne detectors. They are i?j , B^ +7T ^ 2 \ flW, and 

B(x+ 7I 7 2 ) ) corresponding to the "cosine" and "sine" com- 
ponents of the two photocurrents in Eqs. (||) at phases £ 
and x respectively. The modes and B^ are corre- 
lated due to the parametric interaction in Eq. ( pi] ) . This 
interaction, however, does not couple the modes and 

£(A+tt/2)_ 



at the input of the NOPA (i.e., parametric fluorescence). 
In the Schrodinger picture, Eqs. (|ll]) correspond to the 
following state generated at the output of the NOPA: 



|*) = (l-T 2 ) 1/2 5>>''< 



(15) 



where r = tanhr, and the two-mode Fock state |n, m) 
pertaining to Bj and B^ is given by 



B. Photocurrent probability distribution 

Since we are interested in using self-homodyne detec- 
tion to measure the quadratures of two correlated modes 
B^ and , the radio- frequency phase A can always 
be set to zero by shifting the time origin. Then, the 
quadratures X<p(B^) and X^,(B^)) are jointly mea- 
sured. In the following we will use the shorthand nota- 
tion B n = B^ and = X,p(B^), and analogously for 

the input modes = For perfect detectors, the 

joint probability distribution of the "cosine" photocur- 
rents with £ = x m Eqs. (^) coincides with the joint 
probability distribution of the two quadratures X , and 
AT, namely, 



p(x, x'\ </>, ip) = (x, x'; (j>, ip\R\x, x'\ <j>, i/j) 



(12) 



where \x,x';4>,tp) = \x)^ <S> \x')ij> represents the simulta- 
neous eigenvector of the two quadratures X\ and XV 
with eigenvalues x and x' in the Fock space Ht (g> H<-> 

of the two modes -Bj and B^, respectively; and R de- 
notes their joint density operator. For detectors with 
non-unit quantum efficiencies 77 j and , the joint proba- 
bility distribution (%, 4>, ip) of the photocurrents 
is the convolution j240 of the ideal probability in Eq. ([FT 
with Gaussians for each mode of variances 
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(13) 



In this way, the resulting output probability distribution 
can be written in the form 



p njr ,^(x,x';()),ip) = 



1 



TW i?exp 



{x X\f {x> 



2A 2 



2A i 



(14) 



For simplicity, in the following we will assume equal 
quantum efficiencies T]^ — r/^ = r\ for both detectors. No- 
tice that in the limit of unit quantum efficiency, r\ — > 1, 
one has A,, — > 0, and the ideal probability in ([l2]) is 
recovered. 

We are now interested in the simplest case of measure- 
ment, that with u>o ± sidebands in the vacuum state 



H) K 

\Jn\m\ 



(16) 



where |0, 0) denotes the vacuum for both the B \ modes. 
For the state ( |l5| ) the photon-number probability is given 
by 

p(n,m) = \ (n, m\ty)\ 2 

= 5„ m (l-r 2 )r 2n = %f^ T N )" , (17) 
n + 1 \n+ 1 J 

where 



1 - r 2 



v = sinh r 



(18) 



is the average number of photons in each mode at the 
NOPA output due to parametric fluorescence. The main 
feature of the distribution ([[7]) is the perfect correlation 
of the photon numbers in the signal and idler modes. The 
two-mode photon-number probability p(n, m) of Eq. jl^ ) 
is shown in Fig. ^(left). 




FIG. 2. Theoretical two-mode photon-number probability 
p(n,m) of parametric fluorescence in the signal and idler 
(twin-beam) modes (left), and in the ±45°-polarized modes 
(right), given by Eqs. (|Tj) and (^), correspondingly. The 
mean number of photons in each mode n = 10. 

The joint probability distribution of the output pho- 
tocurrents is derived in the Appendix, and is given by 



p n (x,x';<p,ip) 



1 



2ttA 2 



x ( * 



exp 



2A 2 



(x-xir + ( X '-x^r 



*) (19) 
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x exp 



(dl + 4A*)(<P_ 
(x + x') 2 (x 



4A 2 ) 



n'\2 



4A 2 



d 2 +4A 2 d 2 _, 
which can also be cast in the equivalent form 

a K b K 



where 



■ exp 



-a K (x - c K x Y 



K 
j2 



e -i(.*+1>) tanhr 
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i-M 2 ' 



8A 2 



(d 2 +4A 2 )(d 2 
d 2 - d 2 „ 



4A 2 ) 



—re 



6 2 = 



8A 2 
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measure only a single output 
— namely, we ignore the mea- 



In the case that 
photocurrent, say 
sured value of the other photocurrent — the self- 

homodyne detector is equivalent to a conventional ho- 

modyne detector, which measures only the quadrature 

« t 

of mode B^. The output probability distribution is 
given by 



p v {x;cj>) 



2irA 2 



=Tr < pexp 



1^2 



2A 2 



(22) 



implemented to measure the joint distribution of ±45°- 
polarizcd OPA outputs. In that quadrature of the 

annihilation operator 



(20) 



(21) 



B n = iiA^ + vA\ 



(25) 



is detected at a phase arg(a 7r + Ta* ), where the subindex 
7r runs on the two independent ±45° polarizations, 
namely /* and \, and a v is the coherent-state ampli- 
tude of the corresponding central-frequency component 
of the input. Since the interaction ( p5| ) does not couple 
the +45° and —45° modes with each other, the polariza- 
tion non-degenerate OPA is equivalent to two degener- 
ate OPA's. Two-mode joint photon-number distribution 
is just a product of the marginal distributions for each 
mode, and in the case of vacuum-state input sidebands 
is given by @ 



p(n, m) 
p(n, m) 



for n — 2k + 1 or m 
(2fc-l)H (21-1)1! 1 



21 + 1 



2 k + l k\l\ n + 1 

for n — 2k, m = 21 , 



k+l 



(26) 



where the mean photon number h in each mode is 
given by Eq. ( |l8|)_ . The probability distribution ( pfij ) 
is shown in FigTg(right), next to the signal- idler joint 
photon- number distribution of Eq. (jlj) . While the ±45° 
modes exhibit independent photon-number oscillations 
in Fig. H(right), the signal and idler correlations in 
Fig. 0(left) result in oscillations of the total photon num- 
ber. 



III. QUANTUM HOMODYNE TOMOGRAPHY 



where the reduced density operator of the mode Bt is 



= Tr^[|*)<*r] = 



n + 1 \ n + 1 



BIB 



(23) 



with Tr M denoting the partial trace over the Hilbcrt 
space of the undetected mode B^. The reduced den- 
sity operator of the mode Bt in Eq. (^3|) is that of a 
thermal state with the photon-number probability 



p(n) 



n+1 \ n + 1 



(24) 



where the average photon number n is given by Eq. (18) 



The probability distribution of the output photocurrent, 
Eq. (p^), is a Gaussian with variance A 2 = i(n+i) + A 2 , 
centered at zero. This result of self-homodyning of only 
the signal mode has been recently demonstrated experi- 
mentally p3[ . 

Let us note that, while our analysis is aimed at the 
measurement of the joint signal-idler photon distribu- 
tion, a similar self-homodyning approach can also be 



In this section we briefly review the method for re- 
constructing the quantum state that was introduced in 
Refs. [||^] for one field mode. Then we show how it can be 
straightforwardly extended to any number of modes — in 
particular, to the case of two modes involved in the self- 
homodyne detection of the OPA output — and we will ob- 
tain an algorithm similar to those in Refs. Finally, 
we introduce the measurement of the "dressed" state, of- 
ten performed in experiments, as opposed to the "bare" 
state, typically assumed by the theorists. 



A. Single mode detection 

The method for reconstructing the matrix elements of 
the density operator is based on the following resolution 
of the identity on the Hilbert- Schmidt space 



3 = J ^Tr[gD(w)]&(w) , 



(27) 
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where the integral is extended to the complex plane C for 



w. 



and D(w) — exp(— iu*a + wa<) denotes the displace- 
ment operator for the field mode of interest with annihi- 
lation operator a. Equation ( |27| ) simply follows from the 
orthogonality relation for displacement operators 



Tr[D(w)D*(v)] =6 2 (w-v) , 



(28) 



where 62(10) denotes the Dirac delta- function on the com- 
plex plane. By changing to polar variables w = (i/2)ke zc ^, 
Eq. (|27]) becomes 



+ OO 



dk\k\ 



Tr(ge tkx *)e 



(29) 



where Xa = 



i (cre 1 ^ + ae~ 1 ^) denotes the quadrature 
operator for the field mode a. Then we evaluate the 
trace using the eigenvectors {| of and multiply 
and divide the function inside the integral by exp[(l — 
i])k 2 /(8rj)] in the following fashion: 



77 # 



dk\k\ 1-1 1.2 
—e~~ 



/n 7T 

_0 ° 1- 

/ dx P (x,<f>)e lkx e^r k2 e- lkjt * 



(30) 



where p(x,<p) = 4,(x\g\x) $ is the ideal homodyne proba- 
bility. Using the convolution theorem we obtain 



Q = 



dxp v (x; <j))K v (x - X<f,) 



(31) 



where p^(x;^>) is the homodyne probability distribution 
for non-unit quantum efficiency rj, which is the convolu- 
tion of p(x, (f>) with a Gaussian of variance A 2 given in 

Eq. dH). The kernel K v (x - X4,) in Eq. (§|) is formally 
given by 



K v (x - Xa) 
- Re 



dk k exp 



8// 



-k 2 + ik(x - Xa) 



(32) 



where convergence of the integral in Eq. (|32j) for the op- 
erator K^x—Xa) is intended in the weak sense of conver- 
gence of the matrix elements (v\K n {x — Xa)\v') between 
the Hilbert-space vectors \v) and \v'), which are evalu- 
ated before integration. From Eq. (|3l|) it follows that 
the matrix element (v\g\v') can be experimentally ob- 
tained by averaging the function (u|it^(a; — Xa)\v') over 
the quadrature outcomes x that are homodyne detected 
at random phases cj> with respect to the LO, namely 



(v\g\v') = (v^x - Xa)\v') 



(33) 



where the overbar denotes the experimental average. The 
functions {v\K^(x — Xa)\v') for different vectors \v) and 
\v') are called "pattern-functions" after Ref. H. 




FIG. 3. Two-mode photon-number probability p(n, m) of 
the twin-beam state of parametric fluorescence in Eq. (113) 
(two different perspectives), obtained by a Monte-Carlo sim- 
ulation of self-homodyne tomography at unit quantum effi- 
ciency for n = 10 and with 10 6 simulated data. 




FIG. 4. Diagonal elements p(n, n) of Fig. [| (shown by thin 
solid line on an extended abscissa range) with their respective 
error bars in gray shade, compared to the theoretical proba- 
bility (0) (thick solid line) . 

In Ref. 0| the boundness of different types of matrix el- 
ements of the operator kernel K. q (x — Xa) was analyzed 
as a function of the quantum efficiency. It was shown 
that for the photon-number and coherent-state repre- 
sentations these matrix elements become unbounded for 
V < 1/2- The fact that r] — 1/2 is a lower bound for 
measuring the state in any (not exotic) representation 
was thoroughly discussed in Ref. j7j. The way in which 
nonunit quantum efficiency manifests its detrimental ef- 
fect when approaching the lower bound is through in- 
creasingly large statistical errors. Let us restrict our at- 
tention to the photon- number representation. At r\ = 1, 
as proven in Ref. |2j| , the statistical errors of the diago- 
nal matrix elements (n|g|n) saturate at the limiting value 
y/2/N for sufficiently large n, independently of the state 
g (N is the number of data collected in the experiment). 
Also, errors of the off-diagonal elements increase very 
slowly versus the distance from the main diagonal. On 
the other hand, for r/ < 1 the errors increase dramatically 
versus either n or 1 — 77, and eventually become infinite at 
the lower bound r\ = 1/2 p6|. In the next section we will 
see how this behavior manifests itself in the two-mode to- 
mography measurement, on the basis of numerical results 
from Monte-Carlo simulation experiments. 
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B. Multimode detection 



It is easy to see that Eq. (|27j) can be extended because 
of linearity to the case of M modes as follows: 



R 



C M 



lu j2 



M 



R~[[D s (w s ) 



M 



n^K), (34) 



where R now denotes the joint M-mode density opera- 
tor, and £>i(wi) is the displacement operator for the Z-th 
mode. As a consequence, Eq. (33) is extended to the 
multimode measurement in the following way: 



A I 



($\R\$') = (*\]lK m (x l -xf i )\$>) 



(35) 



i=i 



where |$) and |$') are now multimode vectors, and the 
experimental average is taken over the random outcomes 
xi of the joint homodyne measurement of quadratures 
Jtjp , I = l,...,M, of all M modes with random LO 
phases <f>i (we have also let the quantum efficiency to 
be different for each homodyne detector). Equation j35| ) 
agrees with the results obtained in Refs. p8[. 



C. Two-mode tomography through self-homodyning 

As we have seen in Section ||, in the self-homodyne 
measurement one can jointly measure the quadratures 
and X^ of two different modes and _B W , and 
thus, in principle, perform a two-mode tomography of 
the OPA output. However, in order to perform two- mode 
tomography we need uncorrelated phases <f> and ip for the 
quadratures, whereas in the self-homodyne measurement 
they are actually correlated. In fact, one has 



<p = arg(at + ra^) 
ip = arg(a<_> + Tat) 



(36) 



In the case when we are interested in the photon distri- 
bution of one mode only, we can assume — 0. Then, 
by letting the phase of aj fluctuate with a uniform dis- 
tribution from to 2ir, one can perform one-mode to- 
mography of Eq. (|33|). On the other hand, if we are 
interested in the joint photon-number distribution, then 
we can not make the measurement by simply averaging 
the two-mode pattern functions over the experimental 
outcomes as in Eq. (|35|). This is because the LO phases 
<j) and ip in this case are not independent random vari- 
ables. We will show, however, that it is possible to take 
the correlation of <p and ip into account and still per- 
form two-mode tomography by appropriately weighting 
the experimental outcomes in Eq. (|35|). We first rewrite 
Eq. (p5|) in the two-mode case as follows: 



($i = r -± r -± 

K ' ' ' J 2Wo 27T 



dx 



dx' 



xp(x,x' ]( p,ip) {$\K n (x - X^K^x' - X^)\&) . (37) 

We focus our attention on the phase average only. For 
LO's with equal intensities | £*■<-► | = |o!t| and phases tp-t — 
arg(at) and ip^ — arg(a^), one has 



<t> ~ V> = <Pl ~ P-M- , 

(j) + ip = fi + <y5w + 2 arg 



1 + re _i(v t +v " ) 



(38) 



After performing the change of variables 
1 



5 = -x(tpi - <P~) 



(39) 



the average over the phases can be rewritten in terms of 
the average over the sum and difference phases with an 
appropriate weighting function as follows: 



d<p f 27r f +7r dS f 2 * . . fAn . 



where the weighting function is given by 



w(a) 



1 



1 



2vr 1 



2t cos a 



(41) 



Since the input phases tp-r and can certainly be con- 
sidered as random and uncorrelated, the same must hold 
true for their half-sum a and half-difference 8 in Eq. ( |3S|) . 
Then, the measurement of the matrix element in Eq. ( |37| ) 
is obtained by averaging over the experimental random 
phases a and 5 with the weighting function (41). Also, 
the weighting function can be rewritten in terms of the 
gain g(a) of the central-frequency component that is 
given by 



HI 2 



= /i 2 (1 + t 2 + 2tcosct) . (42) 



Hence, the weighting function is simply 

1 



w(a) 



2^7 (<t) 



(43) 



which can be easily and independently measured for ev- 
ery data point while the homodyne data are collected. 

This approach to phase averaging can also be used 
for detec tion of the ±45° modes mentioned in Subsec- 
tion [IB. In that case, the quadrature phases arg(a y + 
Ta^) and arg(a^ + rat ) are independent, but non- 
uniformly distributed. Then, the averaging is done over 
the input phase a — arg(a y) or a = arg(a^), respec- 
tively, with the weighting function ( f43| ) given by the 
phase-sensitive gain of the central component. 
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D. Measuring the "bare" or the "dressed" state 

For non-unit quantum efficiency one can measure the 
density matrix elements for 77 above the bound rj = 1/2. 
However, instead of measuring the density matrix of 
the state R of interest, one can always measure the 
density matrix of the state that has been damped — or 
"dressed" — by the quantum efficiency, without any lim- 
itation for 77, even though such a dressed state would 
be less and less significant for lower quantum efficiencies. 
The concepts of "dressed" and "bare" states are two faces 
of the same measurement description when regarded in 
the equivalent Schrodingcr and Hciscnbcrg pictures. The 
conventional description corresponds to the Heisenberg 
picture, in which the true state — also called the "signal" 
or the "bare" state — is measured and the effect of quan- 
tum efficiency is ascribed to the detector observable (pho- 
tocurrent). In the "dressed" state description, on the 
other hand, one regards the measurement with 77 < 1 on 
the true state R as the corresponding hypothetical "bare 
measurement" with 77 = 1, but now on the "dressed" 
state R n , ascribing the effect of the non-unit quantum 
efficiency to the quantum state itself, rather than to the 
detector. In other words, the effect of the non-unit quan- 
tum efficiency is regarded in a Schrodinger-like picture, 
with the state evolving from R to R^, where the quantum 
efficiency plays the role of a time parameter. 

An easy way to perform tomographic measurement 
on a dressed state is just to use the experimental data 
for 77 < 1 and analyze them using the p atter n func- 
tion with 77 = 1. As shown in Subsection [IB , the ef- 
fect of non-unit quantum efficiency is to convolve the 
quadrature-probability distributions for all LO phases 
with a Gaussian of variance given by Eq. (|l3|). This 
corresponds to convolving the Wigner function with an 
isotropic Gaussian of the same variance in the com- 
plex plane, which, in turn, corresponds to adding Gaus- 
sian noise to the quantum state. In terms of the bare 
state g, the state T n (g) dressed with the Gaussian noise 
is given by [§7l 



f,(e) 



d w 



exp (-\w\ 2 /m) D{w)qD*(w) , (44) 



where the noise-equivalent mean thermal photon number 
m is related to the quantum efficiency through 



m = 2A^ = 



1 — 77 
277 



(45) 



In the multimode case, one needs to apply the transfor- 
mation ([l4|) repeatedly, once per each mode, with the 
corresponding displacement operator of the mode. In 
the context of measuring the bare state, Eq. ( (44| ) was 
exploited in Ref. |28| to show that the measurement is 
possible even in the presence of quantum noise, however, 
with no more than rh — 1/2 thermal photons. 

Another way of dressing the state, which is often em- 
ployed in experimental analysis of the tomographic data 




FIG. 5. Monte-Carlo simulation of self-homodyne tomog- 
raphy of the two-mode photon-number probability p(n, m) for 
quantum efficiency r\ — 0.8 (two different perspectives). The 
state is the same as in Fig . but in its Gaussian-noise dressed 
form Tri(R) [see Eq. (H)]. Notice the smearing effect of the 
non-unit quantum efficiency. 



(see Refs. [y_0|jlljj2^]) is to consider the state that has 
undergone a loss equivalent to 77. In this case, the analy- 
sis is done by rescaling the output photocurrents by ^/rj 
instead of 77 as in Eqs. (||), and then using the pattern 
functions for 77 = 1. ft is easy to see that this procedure 
corresponds to measuring the dressed state A v (g), which 
is related to the bare state p as follows: 



n=0 



(46) 



Again, in the multimode case the transformation ( |46] ) is 
applied separately to all modes. One can also regard the 
state A v (g) in Eq. (|4tj ) as the state of the mode y/r] a + 
y/l — 77 v — instead of the state of just the mode a of 
interest — where v is the independent vacuum-state mode 
responsible for the loss. 

Before concluding this section, we need to say a few 
words regarding the difference between the two dressed 
states A 7) (g) and r i) (g). In the loss model corresponding 

to A v (g), the dressed state loses some signal, and becomes 
the vacuum state in the limit of 77 — ► 0, independently of 
<5, which makes the state less and less meaningful for de- 
creasing 77. On the other hand, in the Gaussian-noise 
model corresponding to T v (q), there is no loss of sig- 
nal, but the state gets an increasingly large number of 
thermal photons for decreasing 77. In this way, the most 
interesting quantum features of the state — as, for exam- 
ple, oscillations in the photon-number probability — are 
lost, as shown in Ref. B, and all states tend to look 
"classical". In the next section we will see these effects 
at work in some Monte-Carlo numerical experiments for 
the two-mode case. 



IV. MONTE-CARLO SIMULATIONS 

In this section we present some numerical results from 
Monte-Carlo simulations of the self-homodyne measure- 
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ment. Our aim is to analyze the feasibility of a real ex- 
periment and to see how many measurements are needed 
for a state reconstruction, especially in presence of the 
detrimental effect of non-unit quantum efficiency of the 
photodetectors. We will restrict our analysis to the mea- 
surement of the joint density matrix R of the two modes, 
B-t and £?<_>, assumed to be in the correlated state given 
by Eq. ©. 

The simulation of the homodyne outcomes is based on 
the probability distribution in Eq. ( ^i"|) , which shows how 
the outcomes can be obtained from a Gaussian random 
generator, starting from the generation of x' , then gener- 
ating x, and finally shifting the latter by c K x' . The phases 
4> and ip of the quadrature are chosen randomly for every 
sample. The density matrix in the photon-number repre- 
sentation is measured by averaging the pattern functions 
over the random data: 

(ni,mi\R\n 2l m 2 ) = 

{ ni \k v (x - X )|n 2 )(mi|£„(a;' - X 4 ,)\m 2 ) . (47) 

The pattern functions for a generic 77 are obtained from 
the pattern functions for rj = 1, using the inverse general- 
ized Bernoulli transformation as in Ref. [^9| . The pattern 
functions for rj = 1, in turn, are obtained from the fac- 
torization formulas of Refs. Q (following our conventions 
for the quadratures, we actually use the factorization for- 
mulas as given in Ref. Q]). In Fig. ||we show the results 
of a simulation for the measurement of the two-mode 
photon-number probability p(n, m) for unit quantum effi- 
ciency. The theoretically expected distribution, given by 
Eq. (|l|), is shown in Fig. g(left). In Fig. f| the diagonal 
elements p(n, n) of Fig. are shown with their respective 
error bars, and compared against the theoretical proba- 
bility of Eq. (0). From both Figs. | and |I] we sec that 
there is an excellent agreement between the theoretically- 
obtained and tomographically-reconstructed joint prob- 
abilities, and the fluctuations in the latter are already 
very small for a number of data samples as low as 10 6 , 
which can be easily acquired within the stability time of 
a typical twin-beam setup. 

In Figs. ^ and ^ the same tomographic measurement 
of Figs. [}] and ^ is reported, but now for a quantum effi- 
ciency 77 = 0.8 for each detector. However, in the recon- 
struction, the pattern functions for 77 — 1 are used. As 
explained in Subsection [III Ej , this corresponds to a mea- 
surement of the state T v (R) that has been dressed by the 
Gaussian-noise equivalent of the quantum efficiency, in- 
stead of a measurement of the true twin-beam state. (For 
values of the quantum efficiency 77 = 0.8 and 77 = 0.9 used 
throughout this paper, the two kinds of state dressing — 
Gaussian-noise or loss — give similar qualitative results.) 
The smearing effect of the non-unit quantum efficiency is 
evident in Fig. |^, where the perfect photon-number cor- 
relation between the two modes is smudged, resulting in 
non-vanishing probabilities p(n, m) for n m. Because 
of the preservation of the normalization in the (n, m) 




n 



FIG. 6. The same as in Fig. but for the diagonal el- 
ements p(n, n) only (thin solid line with error bars in gray 
shade), compared to the theoretical probability ( |l7| ) for the 
bare state (thick solid line). The disagreement between the 
theoretical probability for the bare state R and the simulated 
measurement for the dressed state T V (R) is a typical mani- 
festation of the non-unit quantum efficiency. 




FIG. 7. Reconstruction of the bare state using the pattern 
functions with the correct experimental value of the quan- 
tum efficiency 77 (two different perspectives are shown). Here 
n — 10, rj — 0.9, and we used 10 7 data samples for the 
Monte-Carlo simulation. 

plane, the diagonal p{n, n) is decreased, resulting in the 
evident disagreement in Fig. |^, where the reconstructed 
diagonal elements p(n, n) are reported with relative error 
bars and compared with the theoretical probability ( |T7j ) 
for the bare state. 

In Fig. we present the results of Monte-Carlo simu- 
lation for a realistic measurement of the bare state, but 
now using the pattern functions with the correct exper- 
imental value of the quantum efficiency rj. One can see 
that the smearing effect of the non-unit quantum effi- 
ciency has been cleaned out, which, however, comes at 
the expense of increasing fluctuations for large n. This 
is even more evident in Fig. |^, where the reconstruction 
of the diagonal probability p(n, n) for the bare state is 
shown for two different values, 77 = 0.9 and rj = 0.8, 
of the quantum efficiency. One can see that there is no 
longer the disagreement between the reconstructed and 
the theoretical values, of the kind shown in Fig. g, but 
now the error bars have increased dramatically for larger 
n, becoming worse for smaller 77 [cf. Fig. ^(right)]. 

The off-diagonal number probabilities and the correla- 








n n 

FIG. 8. Reconstruction of the diagonal probability p(n, n) 
for the bare state, using the pattern functions with the correct 
experimental value of the quantum efficiency r\. Here n — 10 
and r\ — 0.9(0.8) in the left (right) figure. The theoretical 
probability (thick solid lines) is superimposed onto the results 
of the Monte-Carlo experiments (10 7 data samples); the latter 
are shown with thin solid lines with statistical errors in gray 
shade. Notice that there is no longer the disagreement shown 
in Fig. ^, but now error bars increase dramatically versus n 
and for smaller r\. 

tion between the two modes can be analyzed by evaluat- 
ing the following sums of matrix elements: 

n 

s(n) = ^2(1, n- l\R\l,n- I) , (48) 

N 

d N (n)= {l,n + l\R\l,n + l) . (49) 

l=max(-n,0) 

The quantity s(n) is the probability distribution for the 
total number of photons in the two modes. The theoret- 
ical result for our state in ( ]lq ) is the oscillating function 

i \ / (1 - t 2 )t 2 ™ n even, , KfV , 
S(n) = \0 nodd, (50) 

similar to the photon-number distribution of a single- 
mode squeezed vacuum j[l],[l7). On the other hand, the 
quantity djy (n) represents the photon- number correlation 
between the two modes, and in the limit N — > oo is the 
Kroneker S n Q for a twin-beam state. For finite N its 
theoretical value for the state in Eq. ( plf ) can be evaluated 
to be 

d N (n) = Ml - r 2 ^ +1 )) . (51) 

In Fig. ^ we show the results of a simulation of the to- 
tal photon- number probability s(n), Eq. (^), for r? = 1 
and compare them to the theoretical value, Eq. (p0|). As 
shown, the theoretically-expected distribution is well re- 
produced from 10 6 data samples with very small sta- 
tistical errors. In Fig. ^ a similar simulation is pre- 
sented as in Fig. ||, but now for a quantum efficiency 
of 7/ = 0.8. The total photon-number probability s(n) 
is reconstructed for both the dressed state and the bare 
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FIG. 9. Oscillations of the total photon-number probabil- 
ity s(n) in Eq. (jl8|) due to the perfect correlation of the pho- 
ton number in the twin-beam state, Eq. @. Thin solid 
line with error bars in gray shade represents the results of a 
Monte-Carlo simulation with unit quantum efficiency, n = 10, 
and 10 6 data samples. Thick solid line is the theoretical re- 
sult, Eq. @. 

state. Once again, one can see the smearing effect of 
the quantum efficiency in the dressed-state case, where 
the oscillations of the total photon number are almost 
completely washed out. On the other hand, the oscil- 
lations are nicely recovered in the reconstruction of the 
bare state, albeit at the expense of increasingly large sta- 
tistical errors. In Fig. hll we present a simulation for 
if = 0.9 to show how these quantum oscillations would 
be detected in an experimentally-feasible measurement 
of the dressed state with n = 4 and 10 6 data samples. 

Regarding measurement of the photon-number cor- 
relation djsr(n) [Eqs. (|4|) and (gj)], comments similar 
to those made for the total photon-number probability 
s(n) hold. Figure [l2] presents the results of a simula- 
tion of the correlation function for the twin-beam state 
with N = n = 10 and unit quantum efficiency, whereas 
Fig. [l^ shows the results of simulations with quantum 
efficiency 7/ = 0.8, once again, reconstructing the correla- 
tion for both the dressed-state and the bare-state cases. 
Here also, the non-unit quantum efficiency in the case of 
dressed-state reconstruction partially smears out the cor- 
relation, which is well recovered in the case of bare-state 
reconstruction. In Fig. E3, we compare the reconstructed 
correlation function for the dressed state in Fig. p^(lcft) 
to that for two modes in uncorrelated coherent states, 
each having the same mean photon number n = 10 as 
the modes of the twin-beam state. One can see that, in 
spite of the detrimental effect of the non-unit quantum 
efficiency, the correlation for the reconstructed dressed 
state is still stronger than that for the uncorrelated coher- 
ent states, the latter representing the standard quantum 
limit. 



V. DISCUSSION 

We have proposed a method for performing two-mode 
optical-homodyne tomography of the twin beams pro- 
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FIG. 10. Similar to Fig. but for a quantum efficiency 
r\ = 0.8. Results for the Gaussian-noise-dressed state recon- 
struction are shown on the left and for the bare state recon- 
struction on the right. Here ii — 10, and we used 5 x 10 6 data 
samples for the left plot and 10 7 data samples for the right 
plot. In the left plot, the ordinate is truncated at the maxi- 
mum value of the simulated probability. The oscillations are 
nicely recovered in the right plot, wherein pattern functions 
with the correct value of quantum efficiency (0.8) were used 
for reconstruction.. 



FIG. 12. Correlation function, Eq. 
state in Eq. (^) with n = 10 and r\ 




FIG. 11. Similar to Fig. 0(left), but for r/ = 0.9, n = 4, 
and 10 6 data samples. 



duced from a nondegenerate optical parametric amplifier. 
The local oscillators (LO's) needed for the homodyne 
tomography are generated in the same parametric pro- 
cess as the twin beams, and, therefore, are automatically 
matched to the signal and idler twin-beam modes. In our 
self-homodyning method, the polarized central spectral 
component at loq serves as the LO for a mode that con- 
sists of two sidebands at ujq ± i7, and the relative optical 
phase between the central component and the sidebands 
can be varied. We have presented a theoretical descrip- 
tion for both one- and two-mode tomography, with main 
focus on measurement of the photon-number distribu- 
tions. For the signal mode alone, a thermal distribution 
of photons is found, in agreement with the results of a 
recent experiment |23| . In the case of two modes, we have 
presented some selected Monte-Carlo simulations of the 
tomographic measurement of the joint photon-number 
distributions, choosing realistic values for the quantum 
efficiency of the photodctcctors. In particular, we have 
analyzed the feasibility of detecting photon-number os- 
cillations and delta-like photon correlation between the 



9|), for the twin beam 
1 reconstructed from 
3 x 10 a data samples. Results of the Monte-Carlo simulation 
(thin solid line with error bars in gray shade) are superim- 
posed onto the theoretical correlation, Eq. (plh, shown by 
thick solid line. 




FIG. 13. Similar to Fig. |l2| but for rj = 0.8. Results for the 
Gaussian-noise-dressed state reconstruction are shown on the 
left and for the bare state reconstruction on the right. In both 
simulations 5 x 10 data samples were used. The non-unit 
quantum efficiency in the dressed-state case partially smears 
out the correlation, which is recovered in the bare-state recon- 
struction (right), however at the expense of increasingly- large 
statistical errors. 



twin-beam modes. We have shown that for ideal pho- 
todetectors such features can be clearly observed even 
with a small number of data samples (10 6 ). However, for 
realistic quantum efficiencies the oscillations are exhib- 
ited with less contrast in the dressed-state reconstruction 
for the same number of data samples. On the other hand, 
for a tomographic measurement of the true output state 
of the OPA, more data samples are needed in order to re- 
duce the statistical errors. Our Monte-Carlo simulations 
show that for a quantum efficiency of rj = 0.9, the oscil- 
lations in the total photon number can be observed, even 
in the dressed-state reconstruction, with as little as 10 6 
data samples, which makes such an experiment feasible. 

We have also shown how the self-homodyning method 
can be used in detection of the ±45°-polarized modes, 
instead of the signal and idler modes. Since in a 
polarization-nondegenerate optical parametric amplifier 
these modes are amplified independently, their joint 
photon-number distribution is factorized into a product 
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FIG. 14. Comparison of the correlation function, Eq. (|l9|), 
for the Gaussian-noise-dressed state of Fig. [Deleft), shown 
by thin solid line with error bars in gray shade, with that for 
two modes in uncorrelated coherent states (thick solid line), 
having the same mean photon number n — 10 per mode. 



of marginal distributions, each exhibiting even-odd oscil- 
lations in its photon number. 

While the focus of our paper has been on the twin- 
beam state, the self-homodyning approach can be ap- 
plied in other instances as well. There are a number of 
mode-matching critical situations where it is possible to 
mix the signal with another mode that underwent a sim- 
ilar generation process. A key requirement in such situa- 
tions would be the scanning of the relative phase between 
the two modes. Among potential applications are detec- 
tion of the superposition (Schrodinger's cat) states, and 
squeezed states generated in optical fibers. 



APPENDIX 



• I, | ,■■ ■ > I ^xx z (x ~\~ X ") z 

(1 - z l ) 1 exp 1 



l-z 2 



(54) 



we can rewrite Eq. (p3|) as 



I*) 



1-K 2 



1/2 



dx 

— oc J —oo 



dx' \x)$ ® \x')^ 



x exp 



Axx'k- {x 2 + x' ){1 + k 2 ) 

~ K 2 



(55) 



where k = r exp[— i(<f> + ip)] (the choice of the branch for 
the square root in the normalization of the state vector 
(|55| ) gives only an overall phase factor that is irrelevant 
for probabilities). Equation ( |55| ) corresponds to the fol- 
lowing joint probability: 



p(x,x';<i>,ip) = 



7t\d K d- 



■ exp 



(x + x') 2 (x-x'f 



(56) 



where d 2 K = |1 + k| 2 /(1 — l K | 2 )- Non-unit quantum ef- 
ficiency of the photodetectors is taken into account by 
evaluating the convolution of the ideal joint probability 
in Eq. (^6|) with Gaussians for each mode of variances 
given by Eq. (|l3|). This immediately leads to Eq. (20). 



ACKNOWLEDGMENTS 

This work was supported in part by the U. S. Office of 
Naval Research. 



In this appendix we derive the joint probability dis- 
tribution, Eq. ( pp[ ) , of the output photocurrents for two- 
mode homodyne detection. 

In the Fock representation, the state at the output of 
the NOPA is given by Eq. dlq), namely, 



|M,) = (l- T 2)V^ T „ Kn) 



(52) 



n=0 



Expanding the Fock state \n) in terms of the quadrature 
representation \x)$ for each mode, one has 



00 r 



00 J — 00 



E 



2™n! 



-H n (V2x)H„(V2x')\x)^ ® \x% , (53) 



where H n (x) denotes the Hermite polynomial of degree 5 
n. Using the following identity p0| ], which is valid for 
any complex number z, 6 



00 / 1 \n 
1 ;Z l 



Y^^-H n {x)H n {x') 



n=0 
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